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ABSTRACT 


Image composition involves the process of embedding a selected region of the 
source image to the target image to produce a new desirable image 
seamlessly. This paper presents an image composition procedure based on 
numerical differentiation using the laplacian operator to obtain the solution of 
the poisson equation. The proposed method employs the red-black strategy to 


speed up the computation by using two acceleration parameters. The method 

is known as modified two-parameter over-relaxation (MTOR) and is an 
Keywords: extension of the existing relaxation methods. The MTOR was extensively 
studied in solving various linear equations, but its usefulness in image 
processing was never explored. Several examples were tested to examine the 
effectiveness of the proposed method in solving the poisson equation for 
MAOR image composition. The results showed that the proposed MTOR performed 
Modified iterative method faster than the existing methods. 
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1. INTRODUCTION 

The image composition process computes the image gradients of both source and target images to 
generate a new desirable image of two composed images seamlessly. Other image processing that employs a 
similar approach includes image matting [1], image stitching [2], surface reconstruction [3], sharpening filter 
[4], image completion [5], and image inpainting [6]. In this paper, we present our work in using numerical 
methods to solve the image composition problem. The Poisson equation is used to model the image 
composition process, in which the image gradients are computed with the Laplacian operator. The problem is 
transformed into a large and often sparse linear system which is suitable to be solved using numerical 
differentiation. In the proposed numerical method, additional acceleration parameters are considered to 
provide more flexibility during parameter tuning and to further improve the execution time. The proposed 
method also employed red-black ordering into the modified variants of the relaxation method, in which pixels 
involved during the composition process are labeled in red or black color nodes. Since red and black nodes 
are computed independently, the computation to obtain the image gradient can be executed in parallel. To 
ensure that all tested methods would generate comparably the same final images, several image similarity 
tests are conducted. 

The use of the Poisson equation in image processing for computing image gradients was introduced 
by Patrick Pérez, Michel Gangnet, and Andrew Blake [7]. Since then several works had been conducted in 
applying image gradients for image processing. R. Raskar, Adrian Ilie, and Jingyi Yu [8] employed an 
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approach based on the gradient domain method which preserves local perceptual properties and minimizes 
usual aliasing, haloing, and ghosting problems. In 2004, Sengupta et al. [1] proposed an approach to the 
natural image matting problem by solving the Poisson equation through the matte gradient domain. It was 
shown in their work that the proposed approach produced encouraging results for various natural images 
which previously were problematic. After that Wang and Yang [2] proposed an image stitching method to 
combines several images that have some overlapped images. A set of the cost function is defined where the 
similarity of the images and visibility of the seam were constructed in the gradient domain. By working in the 
gradient domain rather than the pixel intensities of the images, the authors showed that the proposed 
approach was capable of reducing global inconsistencies and minimizing problems generated by the 
illumination effect. In the work by Gao et al. [5], a completion algorithm based on the gradient domain that 
applied a patch-based technique was reported. This work involved two phases: firstly, they used the gradient 
maps to complete the removed area by applying a patch-based filling algorithm. Secondly, the image was 
then reconstructed from the gradient maps by solving the Poisson equation. To achieve better image 
completion, a matching criterion based on the gradient and color was proposed. 

Image editing based on the gradient domain was also used in the study conducted by Kazhdan et al. 
[3] for surface reconstruction by applying an algebraic approach. In their work, a comprehensive study of 
different integration methods and the problem when working with noisy gradient fields are presented. A 
similar study was proposed by Zhouyu Du, Antonio Robles-Kelly, and Fangfang Lu [9], Dikpal Reddy, Amit 
Agrawal, and Rama Chellappa [10]. Vijayakumar [4] studied the applications of the poisson equation to 
devise several image filters such as sharpening and flattening. They also proposed the solutions in the fourier 
domain which showed a more direct, exact, and efficient approach. This idea was also extended and applied 
to video sequences by employing temporal values as constraints. The image composition approach to reduce 
bleeding artifacts in the generated images was proposed by Michael W. Tao, Micah K. Johnson, and Sylvain 
Paris [11]. The gradient-based method is used for video editing using the variational model by Rida Sadek, 
Gabriele Facciolo, Pablo Arias, et al., [12]. J.-M. Morel, A. B. Petro, and C. Sbert [13] proposed a fourier 
implementation to solve the poisson equation for image processing. Using neumann boundary conditions, the 
poisson equation is solved in the retinex algorithm through fourier implementation by Nicolas Limare, Ana 
Bel en Petro, Catalina Sbert, and Jean-Michel Morel [14]. Zeng et al. [6] proposed a variational framework 
based on a unified four non-local schemes for their work in image inpainting. They applied the gradients 
between patches to generate the inpainted image. The image pyramid method for solving the poisson 
equation was suggested by Hussain and Kamel [15]. Afifi and Hussain [16] proposed the modified Poisson to 
solve the image blending problem. The iterative method to solve the Poisson equation for image editing was 
suggested in [17]. In [18] and [19], the modified over-relaxation and block method were developed. In the 
previous study, the TOR method had been applied to solve path planning problem [20] and many other 
partial differential equations [21]-[24]. 

More recently, Cao et al. [25] proposed a self-embedding image watermarking scheme based on 
reference sharing and the poisson equation. In their work, they established the relationship of each pixel and 
its neighborhood in the original image via Laplacian operator and be converted to compression bits. The 
approach increased the recovered area from the tampered image by reconstructing the relationship between 
each compression bit and each reference bit. In [26], the poisson blending algorithm was applied in the cloud 
removing algorithm for a large area of landsat data. The gradient-domain compositing method is accelerated 
by utilizing a quad-tree approach. The proposed method demonstrated promising results for experimental 
images with cloud coverage of less than 80%. Dongbo Min, Sunghwan Choi, Jiangbo Lu, et al., [27] 
presented an edge-preserving image smoothing method known as fast global smoother. In their method, the 
sparse laplacian matrices are efficiently solved through global objective functions. By combining the efficient 
edge-preserving filters and optimization-based smoothing, the approach obtained comparable runtime to the 
fast edge-preserving filters and also overcame many limitations of the existing local filtering approaches. 
Yang Pan, Mingwu Jin, Shunrong Zhang, et al., [28] proposed an approach that combines deep convolutional 
generative adversarial network (DCGAN) and Poisson blending (PB) to perform image completion tasks. 
Most recently, Nordin Saad, Andang Sunarto, and Azali Saudi [29] employed a modified method by 
combining red-black strategy with accelerated over-relaxation, known as MAOR method, demonstrated 
encouraging results. Works on image editing that involves blending, composition, or completion using the 
Poisson equation are summarized in Table 1. 


Table 1. Works on poisson image editing 


Authors Application Method Authors Application Method 
l ba Successive over- . ; , 
[7] Poisson image editing relaxation (SOR) [19] Image blending SOR and group iteration 
[13] Poisson image editing Fourier transform [26] coe oe ees Gauss-Seidel 
[15] Poisson image editing SOR [28] Map Completion DCGAN and SOR 
[18] Image blending Modified SOR [29] Image composition MAOR 
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2. IMAGE COMPOSITION 

An image is constructed based on a two-dimensional array of pixels as illustrated in Figure 1. In the 
image composition process, the procedure begins by selecting the region of interest R from source image S, 
dR denotes the boundary of R. The selected region is then copied to the desired position of a target image T* 
to produce a new image T. Then, the vector field V of the source image S is computed. Using the Poisson 
equation as a constraint to model the pixels inside the region dR, new pixel values are computed by 
minimizing the gradient difference of V and the selected region of the target image T*. 


The formulation of this minimization problem can be written as: 


ming JJ, IVf — vl? with f Ilao- (1) 


The solution of the minimization problem (1) is the unique solution of the Poisson equation with Dirichlet 
boundary condition. 


In the composition process, the vector field v is a gradient field directly extracted from the source image. 
Therefore, the (2) is rewritten as. 


Af =Vgat0 with flao = f “loo (3) 
where the A denotes the Laplacian operator. 


2.1. Formulation of the proposed method 

To begin the image composition process, the gradient fields in the source image need to be 
computed by solving the Poisson equation that is used to model the problem. Let us consider the standard 
Poisson equation in 2D written as (4). 


ðU . a4U 
axe ta FY) (4) 


The discretization of (4) using the five-point Laplacian operator will result in the following approximation 


(5). 
1 
Ui = (Usi + Ui- j + Ui jy + Uij-1— R fij] (5) 


where h=Ax=Ay. The iterative scheme of (5) is given as (6) [30]. 


(k+1) _ 1 
i,j — 4 


U ju UEO UO USTD — kfl +a- ou (6) 


i+1,j i—1,j i,j+1 i,j—1 i,j 
where œ represents a weighted parameter. In [31], the acceleration over-relaxation (AOR) scheme was 
developed and described as (7). 


(k+1) _ ®777(k) (k) (k) (k) T [rr(k+1) 
UF =a U + Uii, T Ui j4 + Ui j- = h? fi, | + Usa E 
0 


(k+1) (k) (k) 
U EU all eS w)U; j. 


(7) 


where additional parameter, t, is added. An extension to the SOR and AOR is the two-parameter over- 
relaxation (TOR) iterative method [32] which employs two acceleration parameters to provide more 
flexibility and weight tuning choices during the iteration process. The finite-difference approximation 
equation for the TOR method is given as (8). 


(k+1) _ ® frk) (k) (k) (k) T | rr(k+1) 
UF; = [Urki + Uii j + Ui j+ T UF ja g h? fij] + utt 7 


(8) 
(k) Y rrr(k+1) (k) (k) 
Ugra San + (1 — w)U;,, 


The optimal value for the three parameters (œ, t, and y) is in the range between 1 and 2. 
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2.2. Red-black strategy 

Figure 2 illustrates the implementation of a red-black ordering strategy, in which the red and black 
nodes are computed alternatively, thus very suitable for parallel implementation [33]. By applying the 
red-black strategy, the computations of red and black nodes involve two phases. The procedure begins by 
computing red nodes where (1+j) is even. In the second phase, black nodes, where (i+j) is odd, are computed. 
The iteration procedures of these two phases continue until the stopping criterion is satisfied. 
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Figure 1. The standard finite grid of size, N=10 Figure 2. The finite grid with red-black ordering 


Based on the iteration scheme given in (6), the red-black SOR (RBSOR) formula are given as. 


Ue = Slo + Uy + Us + Ue fij] +A- oU for allisjeven (9) 


and 


UD = STUN? + UEP + USED + ULSD — hf] +C- @)U{ for all i+j odd (10) 


The red-black AOR (RBAOR) formula are based on the iteration scheme given in (7) and can be expressed as: 


= male + y™ aie y =a y — h? fi; ;| + (1 — wU for all i+j even (11) 


i+1,j i—1,j i,j+1 i,j—1 
and 
(k+1) _ ®777(k) (k) (k) (k) T 7 7(k+1) (k) 
Uj ag + Uray PU a UO yl ~ Uis t 
(k+1) (k) (k+1) (k) (k+1) (k) (k) f ; 
Visi -7 Ui j + Us j- — Ui j- + Ui 541 7 U il +1- w)U;; for alli + j (12) 
odd 


Lastly, the red-black TOR (RBTOR) formula that is based on (8) are given as. 


(k+1) _ ®777(k) 
UF 4 Lee 


$US) + U, + UO, -Rf yj] + (1 — @) UL? for allitjeven (13) 


i—1,j i,j+1 l,j—1 


and 


ut? = [ue + E E Ue, — bf) += a 


i,J 4 i+1,] i—1,j i,j+1 l,j—1 i—1,j 
(k) (k+1) (k) Y prr(k+1) (k) (k+1) (k) 
Uii + Uit, 7 Uaj] + 4 [Ui j-4 a Ui eh: UF j+ 7” A + (14) 


(L= w)U S for all ¿+j odd 


The computational molecules for red and black nodes are illustrated in Figure 3. Red nodes are calculated 
using the (9), (11), and (13). While, the (10), (12), and (14) are used to calculate black nodes. The 
implementation of these equations is executed iteratively, in which the execution continues until the tolerance 
error rate £, is satisfied. 


Modified iterative method with red-black ordering for image composition using poisson ... (Zakariah Aris) 


20200 ISSN: 2302-9285 





(a) (b) 
Figure 3. Computational molecules for, (a) red, (b) black nodes, respectively 
2.3. Modified iterative methods 


An extension to the RBSOR method is the modified successive over-relaxation (MSOR) [34] that 
can be expressed as: 


I~ E Zu, + U T Ua + Oi =h Es LUA for all i+j even (15) 


and 


Ue) =E jugi + UEP + USIP + UGP — h? fij] + A — BUM for all i+j odd (16) 


i+1,] i—1,j i,j+1 l,j—1 


The parameters a and ß, are fixed such that 0 < a < 2 and 0 < B < 2. The corresponding modified variant of 
RBSOR, namely the modified accelerated over-relaxation (MAOR), are given as. 


y&+D — alee AUV by yO m h? fij] +(1- aU; for alli+jeven (17) 


l,j i+1,j i—1,j i,j+1 i,j—1 


and 


k P yr r(k k k k k 
Ue =E uki t US + Oi toe -kfl + - 


i—1,j i,j+1 i,j—1 i—1,j 
(k) (k+1) (k) (k+1) (k) (k+1) (k) 
Uj + Uit, — Ui j + Ui j- E Ui ja + Ui j+ E eal mak © oe (18) 


B)U;, for all i+j odd 


Lastly, the modified variant of RBTOR, known as the modified two-parameter over-relaxation (MTOR) are 
given as (19) and (20). 
uf) = [yu +uL HUR UR -hft —a@)U® forall itjeven (19) 


i—1,j i,j+1 i,j—1 i,J 


and 


ue? = Fu) + Ue + tu. a 


i,J 4 i+1, i—1,j i,j+1 i,j—1 i—1,j 

(k) (k+1) (k) Y prr(k+1) (k) (k+1) (k) 
Uii j + Uia j B Uaj] + 4 [Ui j-4 T Ui j- + Ui j+ g Uil + (20) 
(1 — p)UÑ for all i+j odd 


Algorithm 1 describes the main procedures that are used in the implementation of the six methods considered 
in this study. 


Algorithm 1: The image composition algorithm for the six methods 


Input : Source image S, Target image G, and pixels inside region dS 
Output : New edited image F 
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/* Initialize U and V (to store the previous and updated pixels inside region dS */ 
initPixels(V, U, S, G) 
/* Initialize parameters a, B, Tt, y */ 
initParams (&, by Ty y) 
while checkConverge(V, U) do 
/* For red nodes (i+j) is even */ 
for i := 1 to n-2 do 
for j := 1 to m-2 do 
if isEven(i+j) then 
V(i,j) := computePixels(U, a, T, y) 
end 
end 
end 


/* For black nodes (i+j) is odd */ 


for i := 1 to n-2 do 
for j := 1 to m-2 do 
if isOdd(itj) then 
V(i,j) := computePixels(U, B, Tt, y) 
end 
end 
end 
end 
F := copyPixels (V) 
return F 


Algorithm 1 is only applied to the pixels involved in the composition process. All other pixels 
outside the region of the composition are ignored during the computation process, thus speed up the overall 
execution time. With MTOR, two acceleration parameters, w, and t, are available for fine-tuning, thus 
provide more choice during the composition process. Additionally, by employing the red-black strategy, two 
additional parameters, a, and B, can be tuned and thus providing additional parameter values option to further 
improve the computational time. 


3. RESULTS AND DISCUSSION 

In this work, three sets of images are used to measure the performance of the considered methods. 
Figure 4 shows all images sets that comprise of source and target images. Table 2 shows the number of pixels 
inside the image masks that were applied for the three image sets. 





SET A (a) SET B (a) 





SET A (b) SET B (b) ” SETC (b) 


Figure 4. These figures are; (a) target, (b) source images for set A, B, and C, respectively 
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Table 2. Number of pixels inside the image masks 


Mask Ballon Goose Eland 
Number of pixels 12.392 7.600 7.939 


3.1. Computational cost measurements 

Table 3, Table 4 and Table 5 show the performance of the six tested methods. The proposed MTOR 
method gives the lowest number of iterations, thus consequently provides the best execution time. All 
modified variants, namely MSOR, MAOR, and MTOR, perform faster than their corresponding standard 
variants of RBSOR, RBAOR, and RBTOR. As shown in the three tables, the modified variants reduce the 
number of iterations by approximately 7% to 14% compared to the standard variants. It is also shown that the 
TOR variants decrease the number of iterations by approximately 8% to 11% and 34% to 43% compared to 
the AOR and SOR variants, respectively. 


Table 3. Computational cost for image set A Table 4. Computational cost for image set B 
METHODS Number of iterations Execution time METHODS Number of iterations Execution time 
RBSOR 1404 19s 271ms RBSOR 488 5s 893ms 
RBAOR 1134 15s 507ms RBAOR 389 4s 748ms 
RBTOR 1046 14s 342ms RBTOR 357 4s 393ms 
MSOR 1298 18s 759ms MSOR 450 5s 493ms 
MAOR 1019 16s 189ms MAOR 348 4s 240ms 
MTOR 928 13s 70ms MTOR 315 3s 818ms 


Table 5. Computational cost for image set C 


METHODS Number of iterations Execution time 


RBSOR 575 6s 977ms 
RBAOR 465 5s 674ms 
RBTOR 429 5s 442ms 
MSOR 531 7s 457ms 
MAOR 418 5s 866ms 
MTOR 380 5s 355ms 


It is clearly shown from the results shown in Table 3, Table 4 and Table 5 that the MTOR 
outperforms the other methods. All variants of TOR are superior to their AOR and SOR variants, in which 
the modified variants performed slightly faster than their corresponding standard variants. Also, since all 
methods applied red-black ordering, the proposed methods can be implemented for parallel execution. 


3.2. Image quality measurements 

The image quality measurement process is used to compare the similarity between the final output 
images. Based on the statistical method analysis of variance (ANOVA) [35], seven metrics are used namely 
mean square error (MSE), structural similarity index (SSIM) [36], structural content (SC), normalized cross- 
correlation (NCC), normalized absolute error (NAE), average difference (AD) and maximum difference 
(MD) [37]. The simple MSE can be written as (19). 


1 
where A and B represent the pixel values of reference and target images, and (m x n) is the size of the image. 
MSE is used to measures the difference between pixel values in A and B, in which a smaller value means 
higher similarity. The ideal MSE value of 0 is obtained when the two images A and B are identical. Another 
similarity test that can be used to compare the two images is SSIM that is written as (20) [35]. 


(2PxPytCyz)(2TxytC2) 
(PEt 95 +C1)(TL+TH+C2) 


SSIM(x,y) = (20) 


where @x, Ọy and T, ty denote the mean intensity and standard deviation set of image block x and image block 
y, respectively, while Txy denote their cross-correlation. Cı and C2 are small constants value to avoid 
instability problems when the denominator is too close to zero [38]. If the obtained SSIM is equal to 1, it 
means the two images identical. Similarly, SC and NCC values that are equal to 1 can be used to indicate that 
the two images are identical. SC value can be obtained using the (21). 
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NCC measurement compares the two images and is given as (22). 
Ajj XBij 
NCC = Vita Divi (22) 
ij 


In contrast, NAE and AD values that are close to 0 means the two images are identical. Both values can be 
obtained as. 


m n 
Vint X j=1(l4ij-Bij 


NAE = 
A Eji)? 


(23) 


and 
1 
AD = — dizi Lj=1lAij — Bij], (24) 


where A and B are the reference and target images, respectively. Lastly, the MD value that is used to obtain 
the maximum difference of pixel value between the two images A and B. MD 1s given the (25). 


If the obtained MD value is very small, it means the similarity between the two images 1s higher. All 
these seven metrics are used to compare the similarity of the output images generated by the six tested 
methods. As shown in Figure 5, Figure 6, and Figure 7, the generated images are all visually identical. The 
results of the seven metrics measurement are given in Table 6, Table 7 and Table 8. 


| | = 


$; 


. ’ 


= tS | 












Figure 5. The output images of set A 
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Figure 6. The output images of set B 





Figure 7. The output images of set C 


Based on the similarity measurement shown in Table 6, Table 7 and Table 8, the MSE, NAE, and 
AD values are approximately near to 0. Meanwhile, the SSIM and NCC are very close to 1, which means the 
generated images are all almost identical. The maximum difference (MD) in pixel values is very small 
(between 2 to 8). The obtained values show that the six methods produce identical images. 

The results obtained by the considered algorithms show that red-black strategy provides an 
alternative to the existing approaches for solving image composition problem. By applying a chessboard like 
point coloring, the red-black strategy is very suitable for parallel implementation, since the computation for 
red and black nodes can be done independently. 


Table 6. Similarity measurement for image set A 
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METHODS SSIM MSE SC NCC NAE AD MD 
RBSOR 0.99993 0.23408 0.99912 1.00044 0.00039 0.059655 7 
RBAOR 0.99992 0.27711 0.99904 1.00047 0.00042 0.064970 8 
RBTOR 0.9999 1 0.28855 0.99902 1.00048 0.00043 0.0663 10 8 

MSOR 0.99992 0.24488 0.99910 1.00045 0.00039 0.060975 a 
MAOR 0.99991 0.28701 0.99902 1.00048 0.00043 0.066140 8 
MTOR 0.9999 1 0.29958 0.99900 1.00049 0.00044 0.067555 8 
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Table 7. Similarity measurement for image set B 


METHODS SSIM MSE SC NCC NAE AD MD 
RBSOR 0.99999 0.00591 0.99995 1.00002 0.00003 0.005580 2 
RBAOR 0.99999 0.00646 0.99995 1.00003 0.00003 0.006050 3 
RBTOR 0.99999 0.00667 0.99995 1.00003 0.00003 0.006235 3 
MSOR 0.99999 0.00601 0.99995 1.00002 0.00003 0.005670 2 
MAOR 0.99999 0.00665 0.99995 1.00003 0.00003 0.006220 3 
MTOR 0.99999 0.00678 0.99994 1.00003 0.00003 0.006330 3 


Table 8. Similarity measurement for image set C 


4. 


METHODS SSIM MSE SC NCC NAE AD MD 
RBSOR 0.99999 0.11480 0.99900 1.00050 0.00039 0.049680 4 
RBAOR 0.99998 0.13669 0.99890 1.00055 0.00042 0.054495 4 
RBTOR 0.99998 0.14307 0.99888 1.00056 0.00043 0.055795 5 

MSOR 0.99999 0.11947 0.99898 1.00051 0.00039 0.050710 4 
MAOR 0.99998 0.14243 0.99888 1.00056 0.00043 0.055670 5 
MTOR 0.99998 0.14919 0.99885 1.00057 0.00044 0.057025 5 

CONCLUSION 


For solving the image composition problem, the red-black ordering implementation is used for the 


six iterative methods considered in this study. The experimental results show that the MTOR method 
provides the fastest rate of convergence compared to the previous MSOR, MAOR, and all Red-Black 
variants. Most of the previous works did not utilize the red-black strategy that supports parallelism. Seven 
image similarity indexes were used to compare the generated images. The similarity tests conducted on the 
images show that the suggested methods obtained almost identical images with no noticeable differences. In 
the future, the work can be extended by applying more advanced point iterative methods such as half- and 
quarter-sweep iteration approaches. The iterative methods based on block can also be considered. 
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